Template for GeoPressureR: Untitled

Show code

library(GeoPressureR)
library(tidyverse)
library(leaflet)
library(leaflet.extras)
library(raster)
library(dplyr)
library(ggplot2)
library(plotly)
knitr::opts_chunk$set(echo = FALSE)
load(paste0("../data/1_pressure//", params$gdl_id, "_pressure_prob.Rdata"))
load(paste0("../data/2_light/", params$gdl_id, "_light_prob.Rdata"))
load(paste0("../data/3_static/", params$gdl_id, "_static_prob.Rdata"))
# load(paste0("../data/4_basic_graph/", params$gdl_id, "_basic_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_wind_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_grl.Rdata"))
col <- rep(RColorBrewer::brewer.pal(8, "Dark2"), times = ceiling(max(tag$sta$sta_id) / 8))

Altitude

Altitudes are computed based on pressure measurement of the geolocation, corrected based on the assumed location of the shortest path. This correction accounts therefore for the natural variation of pressure as estimated by ERA-5. The vertical lines indicate the sunrise (dashed) and sunset (solid).

Show code

p <- ggplot() +
  # geom_line(data = tag$pressure, aes(x = date, y = obs), colour = "grey") +
  geom_line(data = do.call("rbind", shortest_path_timeserie), aes(x = date, y = altitude)) +
  geom_line(data = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0), aes(x = date, y = altitude, col = factor(sta_id))) +
  geom_vline(data = twl, aes(xintercept = twilight, linetype = ifelse(rise, "dashed", "solid"), color="grey"), lwd=0.1) +
  theme_bw() +
  scale_colour_manual(values = col) +
  scale_y_continuous(name = "Altitude (m)")

ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)

Wintering location

Show code

file <- paste0("figure_print/wintering_location/wintering_location_",params$gdl_id,".png")
if(file.exists(file)){
  knitr::include_graphics(file)
}

Latitude time

Show code

 tmp <- lapply(pressure_prob, function(x) {
    mt <- metadata(x)
    df <- data.frame(
      start = mt$temporal_extent[1],
      end = mt$temporal_extent[2],
      sta_id = mt$sta_id
    )
  })
  tmp2 <- do.call("rbind", tmp)

sim_lat <- as.data.frame(t(path_sim$lat)) %>%
  mutate(sta_id = path_sim$sta_id) %>%
  pivot_longer(-c(sta_id)) %>%
  left_join(tmp2,by="sta_id")

sim_lat_p <- sim_lat %>%
  filter(sta_id==max(sta_id)) %>%
  mutate(start=end) %>%
  rbind(sim_lat)

sp_lat <- as.data.frame(shortest_path) %>% left_join(tmp2,by="sta_id")

sp_lat_p <- sp_lat %>%
  filter(sta_id==max(sta_id)) %>%
  mutate(start=end) %>%
  rbind(sp_lat)

p <- ggplot() +
  geom_step(data=sim_lat_p, aes(x=start, y=value, group=name), alpha=.07) +
  geom_point(data=sp_lat_p, aes(x=start, y=lat)) +
  xlab('Date') +
  ylab('Latitude') +
  theme_light()

ggplotly(p, dynamicTicks = T)

Shortest path and simulated path

The large circles indicates the shortest path (overall most likely trajectory) estimated by the graph approach. The size is proportional to the duration of stay. The small dots and grey lines represents 10 possible trajeectories of the bird according to the model.

Click on the full-screen mode button on the top-left of the map to see more details on the map.

Show code

sta_duration <- unlist(lapply(static_prob_marginal, function(x) {
  as.numeric(difftime(metadata(x)$temporal_extent[2], metadata(x)$temporal_extent[1], units = "days"))
}))
pal <- colorFactor(col, as.factor(seq_len(length(col))))
m <- leaflet(width = "100%") %>%
  addProviderTiles(providers$Stamen.TerrainBackground) %>%
  addFullscreenControl() %>%
  addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
  addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = pal(factor(shortest_path$sta_id, levels = tag$sta$sta_id)), weight = sta_duration^(0.3) * 10)

for (i in seq_len(nrow(path_sim$lon))) {
  m <- m %>%
    addPolylines(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = 0.5, weight = 1, color = "#808080") %>%
    addCircles(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = .7, weight = 1, color = pal(factor(shortest_path$sta_id, levels = tag$sta$sta_id)))
}
m

Marginal probability map

The marginal probability map estimate the overall probability of position at each stationary period regardless of the trajectory taken by the bird. It is the most useful quantification of the uncertainty of the position of the bird.

Show code

li_s <- list()
l <- leaflet(width = "100%") %>%
  addProviderTiles(providers$Stamen.TerrainBackground) %>%
  addFullscreenControl()
for (i_r in seq_len(length(static_prob_marginal))) {
  i_s <- metadata(static_prob_marginal[[i_r]])$sta_id
  info <- metadata(static_prob_marginal[[i_r]])$temporal_extent
  info_str <- paste0(i_s, " | ", info[1], "->", info[2])
  li_s <- append(li_s, info_str)
  l <- l %>%
    addRasterImage(static_prob_marginal[[i_r]], colors = "OrRd", opacity = 0.8, group = info_str) %>%
    addCircles(lng = shortest_path$lon[i_s], lat = shortest_path$lat[i_s], opacity = 1, color = "#000", weight = 10, group = info_str)
}
l %>%
  addLayersControl(
    overlayGroups = li_s,
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  hideGroup(tail(li_s, length(li_s) - 1))

Wind assistance

Show code

  fun_marker_color <- function(norm){
    if (norm < 20){
      "darkpurple"
    } else if (norm < 35){
      "darkblue"
    } else if (norm < 50){
      "lightblue"
    } else if (norm < 60){
      "lightgreen"
    } else if (norm < 80){
      "yellow"
    } else if (norm < 100){
      "lightred"
    } else {
      "darkred"
    }
  }
  fun_NSEW <- function(angle){
    angle <- angle  %% (pi* 2)
    angle <- angle*180/pi
    if (angle < 45/2){
      "E"
    } else if (angle < 45*3/2){
      "NE"
    } else if (angle < 45*5/2){
      "N"
    } else if (angle < 45*7/2){
      "NW"
    } else if (angle < 45*9/2){
      "W"
    } else if (angle < 45*11/2){
      "SW"
    } else if (angle < 45*13/2){
      "S"
    }else if (angle < 45*15/2){
      "SE"
    } else {
      "E"
    }
  }

  sta_duration <- unlist(lapply(static_prob_marginal,function(x){as.numeric(difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1],units="days"))}))

  m <-leaflet(width = "100%") %>%
    addProviderTiles(providers$Stamen.TerrainBackground) %>%  addFullscreenControl() %>%
    addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
    addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#000", weight = sta_duration^(0.3)*10)

  for (i_s in seq_len(grl$sz[3]-1)){
    if (grl$flight_duration[i_s]>5){
      edge <- which(grl$s == shortest_path$id[i_s] & grl$t == shortest_path$id[i_s+1])

      label = paste0( i_s,': ', grl$flight[[i_s]]$start, " - ", grl$flight[[i_s]]$end, "<br>",
                      "F. dur.: ", round(grl$flight_duration[i_s]), ' h <br>',
                      "GS: ", round(abs(grl$gs[edge])), ' km/h, ',fun_NSEW(Arg(grl$gs[edge])),'<br>',
                      "WS: ", round(abs(grl$ws[edge])), ' km/h, ',fun_NSEW(Arg(grl$ws[edge])),'<br>',
                      "AS: ", round(abs(grl$as[edge])), ' km/h, ',fun_NSEW(Arg(grl$as[edge])),'<br>')

      iconArrow <- makeAwesomeIcon(icon = "arrow-up",
                                   library = "fa",
                                   iconColor = "#FFF",
                                   iconRotate = (90 - Arg(grl$ws[edge])/pi*180) %% 360,
                                   squareMarker = TRUE,
                                   markerColor = fun_marker_color(abs(grl$ws[edge])))

      m <- m %>% addAwesomeMarkers(lng = (shortest_path$lon[i_s] + shortest_path$lon[i_s+1])/2,
                                   lat = (shortest_path$lat[i_s] + shortest_path$lat[i_s+1])/2,
                                   icon = iconArrow, popup = label)
    }
  }
  m

Histogram of Speed

Show code

edge <- t(graph_path2edge(path_sim$id, grl))
nj <- ncol(edge)
nsta <- ncol(path_sim$lon)

speed_df <- data.frame(
  as = abs(grl$as[edge]),
  gs = abs(grl$gs[edge]),
  ws = abs(grl$ws[edge]),
  sta_id_s = rep(head(grl$sta_id,-1), nj),
  sta_id_t = rep(tail(grl$sta_id,-1), nj),
  flight_duration = rep(head(grl$flight_duration,-1), nj),
  dist = geosphere::distGeo(
    cbind(as.vector(t(path_sim$lon[,1:nsta-1])), as.vector(t(path_sim$lat[,1:nsta-1]))),
    cbind(as.vector(t(path_sim$lon[,2:nsta])),   as.vector(t(path_sim$lat[,2:nsta])))
  ) / 1000
) %>% mutate(
  name = paste(sta_id_s,sta_id_t, sep="-")
)

plot1 <- ggplot(speed_df, aes(reorder(name, sta_id_s), gs)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot2 <- ggplot(speed_df, aes(reorder(name, sta_id_s), ws)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot3 <- ggplot(speed_df, aes(reorder(name, sta_id_s), as)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot4 <- ggplot(speed_df, aes(reorder(name, sta_id_s), flight_duration)) + geom_point() + theme_bw() +scale_x_discrete(name = "")

subplot(ggplotly(plot1), ggplotly(plot2), ggplotly(plot3), ggplotly(plot4), nrows=4, titleY=T)

Table of transition

Show code

alt_df = do.call("rbind", shortest_path_timeserie) %>%
    arrange(date) %>%
    mutate(
      sta_id_s = cummax(sta_id),
      sta_id_t = sta_id_s+1
    ) %>%
    filter(sta_id == 0 & sta_id_s > 0 ) %>%
    group_by(sta_id_s, sta_id_t) %>%
    summarise(
      alt_min = min(altitude),
      alt_max = max(altitude),
      alt_mean = mean(altitude),
      alt_med = median(altitude),
    )

  trans_df <- speed_df  %>%
    group_by(sta_id_s,sta_id_t,flight_duration) %>%
    summarise(
      as_m = mean(as),
      as_s = sd(as),
      gs_m = mean(gs),
      gs_s = sd(gs),
      ws_m = mean(ws),
      ws_s = sd(ws),
      dist_m = mean(dist),
      dist_s = sd(dist)
    ) %>%
    left_join(alt_df)

trans_df %>% kable()

sta_id_s	sta_id_t	flight_duration	as_m	as_s	gs_m	gs_s	ws_m	ws_s	dist_m	dist_s	alt_min	alt_max	alt_mean	alt_med
1	2	3.4166667	34.57011	17.494468	47.68620	21.711983	18.003929	3.8248656	162.92785	74.18261	157.25331	1060.0105	654.1848	604.5139
2	3	7.3333333	37.18249	15.983395	44.39281	16.646386	11.408376	5.5214162	325.54730	122.07350	240.07618	985.7225	595.2773	547.5521
3	4	7.9166667	40.58991	15.977018	56.06959	18.303820	19.826997	5.5771796	443.88426	144.90524	-230.29992	1404.8480	298.0325	214.8520
4	5	4.9166667	33.81206	17.370380	34.80925	16.913821	4.123158	2.5650512	171.14546	83.15962	22.06228	1241.4214	487.8901	506.4481
5	6	9.8333333	41.63452	11.463145	50.63115	16.145250	11.780202	2.6895206	497.87295	158.76162	9.47344	1670.3913	679.8558	587.9694
6	7	9.4166667	36.63912	13.033330	36.98587	12.525579	12.530809	5.0572737	348.28359	117.94921	127.30825	3458.0496	1838.6165	1684.6414
7	9	28.4166667	55.26767	17.586456	64.81036	14.110518	11.531872	3.7032344	1841.69430	400.97388	NA	NA	NA	NA
9	10	8.5000000	62.58095	19.842185	71.09937	18.510293	16.076736	2.4089642	604.34461	157.33749	935.15268	2205.3068	1804.8507	1854.0280
10	11	8.9166667	41.04600	17.261752	48.03535	18.698191	20.129450	4.1477964	428.31520	166.72553	658.94857	3100.9877	1691.4406	1707.8348
11	12	5.2500000	41.59498	21.239332	40.42210	19.901930	3.833050	1.4844871	212.21603	104.48513	295.45738	1527.9584	1134.8679	1336.0644
12	13	7.1666667	44.52462	17.710289	38.32131	17.694810	8.304144	1.8732288	274.63602	126.81281	359.33122	1215.7742	706.7753	650.7422
13	14	0.8333333	35.96519	16.702519	36.62206	20.024884	9.981134	3.0560772	30.51839	16.68740	136.05294	219.7510	177.9019	177.9019
14	16	12.6666667	12.42522	6.712675	12.08427	7.529162	3.545367	0.9657721	153.06746	95.36939	NA	NA	NA	NA
16	17	3.7500000	32.95460	11.883280	32.75673	11.478415	2.103072	1.1820068	122.83775	43.04406	220.08353	411.3244	317.3691	321.5518
17	18	1.5000000	33.64262	10.438046	39.00593	17.020175	17.288656	1.6683386	58.50889	25.53026	141.33298	171.3591	154.2831	150.1572
18	19	10.6666667	37.47178	11.796753	43.69214	10.826953	18.739869	2.9062068	466.04952	115.48749	243.59053	1672.0438	896.6158	895.9195
19	20	4.9166667	33.69398	14.300953	42.22981	17.339492	19.528778	4.9322794	207.62988	85.25250	312.52187	1782.4563	906.9879	825.7417
20	21	4.2500000	36.99632	14.036391	28.68552	18.504806	25.468923	4.2639659	121.91347	78.64543	287.69319	1124.9316	655.0057	669.5111
21	22	10.6666667	36.38271	17.266830	30.44620	12.582570	16.071009	7.0997736	324.75946	134.21408	432.44483	1905.1577	794.3203	628.2020
22	23	10.5000000	39.79834	16.936986	38.94978	16.247836	12.807723	5.7257082	408.97265	170.60228	789.05124	3154.8776	2129.4470	1842.2562
23	24	10.2500000	39.65864	13.748037	78.32582	15.589634	43.528736	5.2998547	802.83962	159.79375	1035.30185	3786.9226	3129.5720	3261.6373
24	26	23.5000000	51.53299	8.173692	113.85681	3.927738	66.987229	6.5230024	2675.63511	92.30185	NA	NA	NA	NA
26	27	7.9166667	26.75544	7.335790	31.91815	8.167588	8.769154	1.2363085	252.68532	64.66007	340.90962	1229.2793	706.1691	737.0251
27	28	5.7500000	39.82505	16.393176	57.14521	21.375869	20.570207	6.4380061	328.58494	122.91125	474.38752	2804.7508	1480.3557	1389.4821
28	29	6.5000000	61.13500	15.352408	62.17359	17.253713	25.230900	5.6896699	404.12837	112.14913	207.11052	2572.7330	964.6875	450.2296

18LX

Altitude

Wintering location

Latitude time

Shortest path and simulated path

Marginal probability map

Wind assistance

Histogram of Speed

Table of transition